Detection of confinement and jumps in single molecule membrane trajectories 
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We propose a novel variant of the algorithm by Simson et al. [R. Simson, et al, Biophys. J. 69, 
989 (1995)]. Their algorithm was developed to detect transient confinement zones in experimental 
single particle tracking trajectories of diffusing membrane proteins or lipids. We show that our 
algorithm is able to detect confinement in a wider class of confining potential shapes than Simson et 
al.'s one. Furthermore it enables to detect not only temporary confinement but also jumps between 
confinement zones. Jumps are predicted by membrane skeleton fence and picket models. In the case 
of experimental trajectories of /^-opioid receptors, which belong to the family of G-protein-coupled 
receptors involved in a signal transduction pathway, this algorithm confirms that confinement cannot 
be explained solely by rigid fences. 



One of the central issues of contemporary cellular bi- 
ology is to establish the relationship between dynamical 
organization and biological functions of membrane con- 
stituents. The development of single particle tracking 
(SPT) techniques gives enthusiastic new insights into dy- 
namical organization of membranes, that were so far in- 
accessible by ensemble-average methods. Indeed, the dif- 
fusive motion of molecules of interest (proteins or lipids) 
at the surface of living cells can be followed with a nano- 
metric resolution, after labeling them by means of fluo- 
rophores, gold colloids, latex beads, or quantum dots 0. 
Yet, new specific and performant tools must be developed 
to extract valuable information from these trajectories. 

After more than 15 years of efforts using SPT tech- 
niques, the question of the cell membrane organization 
and compartmentalization is still the matter of intense 
and controversial debate H H S S H 0- What is 
the origin of the confinement quite generally observed 
in tracking experiments? How is it related to the trans- 
mission of signals through the membrane? Confinement 
is indeed commonly observed in SPT trajectories: the 
diffusive motion is not purely Brownian, but rather af- 
fected by either rigid obstacles, confining domains such 
as rafts (or other signaling platforms), or molecular inter- 
actions. Confinement can be transient Q, the molecule 
being now and then trapped in "transient confinement 
zones" (TCZ). Experimental situations also exist where 
the molecule is always confined, while showing a long- 
term slow diffusion. There exist different models that 
account for such a behavior. In the "membrane skeleton 
fence and picket models" 0, |^ , the confinement is due to 
the cytoskeleton of the cell close to the membrane, or by 
proteins anchored to it, which form rigid corrals. Succes- 
sive hops between adjacent domains result in a slow long- 
term diffusion of the molecules (Fig. The recent al- 
ternative "interacting Brownian particle model" pro- 
poses that - in the case studied in this reference - barri- 
ers do not satisfyingly explain the observed confinement. 



which more likely originates from long-range attractive 
interactions between membrane proteins. The latter form 
auto-assembled aggregates in which proteins are trapped. 
The long-term diffusive behavior is the manifestation of 
the diffusion of the center of mass of the assembly. We 
demonstrate that the algorithm studied in this paper is 
able to discriminate between these models. Indeed, in 
a fence and picket model, proteins regularly jump from 
a confining zone to an adjacent one. By contrast, the 
"interacting Brownian particle model" does not require 
jumps to account for long-term slow diffusion. We show 
that, despite statistical fluctuations, our algorithm de- 
tects jumps with good confidence, when they exist. We 
calculate how many jumps are theoretically expected in 
a fence and picket model, and we compare this number 
to the effectively detected ones in experimental trajecto- 
ries of Ref. 0] (/Lt-opioid receptors). We find that there is 
an unequivocal discrepancy. This proves that rigid fences 
cannot be considered as the unique source of confinement. 

Beyond this particular example, our algorithm intends 
to be applicable to a wide range of experimental situa- 
tions. It responds to an increasing demand consecutive to 
the rapid development of SPT experiments. It is a chal- 
lenge to develop a simple and reliable tool to discriminate 
between different sources of confinement or more sim- 
ply between confined and non-confined trajectories. The 
present Letter intends to propose such a robust tool. 

Simson et aL's algorithm [l^l has been designed to de- 
tect transient confinement. It is based on the following 
principle. Consider a Brownian trajectory r(t) on a time 
interval [to, to + St]. The maximum of ||r(t) — r(to)|p 
on [to, to -I- St], denoted by r^axj scales like DSt where 
D is the diffusion constant. Then the authors define a 
confinement index A (denoted by L in their paper) that 
is an affine function of DSt/r^^,^. D is determined by 
measuring the slope at the origin of the mean square de- 
viation MSD(t). If the diffusion is confined in a domain 



of typical size L, then 



is limited by L , and A is 
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larger than in the case of free Brownian diffusion. The 
authors determine a threshold Ac- Roughly speaking, if 
A > Xr, the diffusion is confined, otherwise it is free (see 
Ref. im for more details). Along a trajectory, X{t) is 
calculated on sliding intervals [t — dt/2,t + St/2]. The 
plot of X{t) indicates TCZs as intervals where A > Ac. 

We show that even though Simson et al's method is 
applicable to a large variety of experimental cases, there 
exist situations of great interest where it is not opera- 
tional. First, artifactual detections of TCZs can happen 
when D varies along the trajectory This problem is 
fixed by computing D locally by the same method, on in- 
tervals of duration a few seconds. More importantly, this 
method fails in detecting confinement in non-flat poten- 
tials. For example, in a quadratic well of typical width L 
at temperature T, the molecule is likely to explore regions 
of energy of several ksT where r ^ L, and the measured 
r^j^x fluctuates a lot around its typical value, depending 
on whether such rare points are in the trajectory or not. 
We experienced that it happens that Simson et al.'s algo- 
rithm does not detect flagrant confinement in quadratic 
potentials or in experimental trajectories (see Fig. [3 for 
an example). The algorithm proposed in Ref. jllj, also 
based on r^^^^, presents the same limitations. This is 
intrinsic to the methods, namely the choice of r^^^^ to 
characterize trajectory wanders, and not to a particular 
choice of parameters. 

For this reason, we modify the above algorithm as fol- 
lows: instead of calculating r^^^, we compute the vari- 
ance Ar2(t) of r on intervals [t - 5t/2,t + St/2]. Rare 
points wandering far away from the potential minimum 
thus have a low weight in Ar^ (t) , which gives a more ac- 
curate measurement of the typical width of the confining 
potential. Of course, in the case of flat potentials de- 
limited by rigid fences, both algorithms present a similar 
efficiency. A novel confinement index is now defined 



Up to a numerical prefactor, A is the ratio of the vari- 
ance of a free random walk to the one of the trajec- 
tory under study. If it is unconfined, A will be of or- 
der unity, whereas it will be large in the converse case. 
We calculate the typical values of A in the respective 
cases of free Brownian two-dimensional trajectories and 
confined ones. We model our Brownian molecule by an 
over-damped Langevin particle: dr/dt = rjit), where t] 
is a Gaussian white noise: (rj) = and {r]i{t)'>lj {t')) — 
2DSijS{t ~ t'). Then the mean position r and the mean 
square position are calculated on a single trajectory 
before being averaged over noise, leading to Ar^ = (r^ — 
r^) — {2/3)DSt. Statistical fluctuations of this quantity 
can also be calculated using Wick's theorem |0| to com- 
pute 4-time correlators of r;: A(r2 — r^) = j^DSt. Thus 
Asrotim ^ 1/(2/3 — 2\/2/3\/5) — 4 for a pure Brownian 



trajectory, independently of St and D, as checked on nu- 
merical trajectories. Now we consider confined diffusion 
in a square box of side L. If one averages over N 1 
independent positions, one gets Ar^ = /6, 

Aeon/ = 6DSt/L^. (2) 

In this confined case, the statistical fluctuations of Ar^ 
vanish at large N. Hence TCZs will be distinguishable 
from pure Brownian trajectories if Aconf 3> ^Brown i-e- 
if Aconf > 4 or DSt/L^ > 2/3. As expected, St must 
be sufficiently large to enable detection. Note that above 
we have asked for the number of independent images to 
be sufficiently large. In practice, positions in successive 
images are correlated, because the equilibriation time to 
explore a box of side L is r = L^/n^D H|. One must 
make sure that ^ I with respect to this time. 
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FIG. 1: (Color online) Example of trajectory in a square grid 
of rigid barriers of periodicity L. At short times, the particle 
diffuses in a closed box, where it stays on average a time Tres • 
Periodically, it jumps from one box to an adjacent one, thus 
resulting in a slow long-term diffusion with constant Dm- 

Our new algorithm detects confinement with excellent 
reliability, whatever the shape of the box or of any non- 
flat confining potential. Now how does it detects jumps 
over fences? We suppose first that the molecule evolves 
in a grid of square boxes of side L, separated by rigid 
barriers, and that it can jump over a barrier by thermal 
activation (Fig.^. It has a long-term slow diffusion with 
a constant Dm- Ar^ — ADm t at large t. If r^es is the 
average time of residence in a box, i. e. the average time 
between successive jumps, then = ADMTres- When 
the molecule is confined in a single box, A = Aeon/- If 
there is a jump at time io in a segment /= [t — St/2,t + 
St/ 2], the particle is virtually in a larger box 2L x L: it 
spends a time tQ—t+St/2 in a box and the remaining time 
in an adjacent one. liSt ^ r, the probability distribution 
of the molecule position is close to uniform in each box. 
This permits to calculate Ar^ by partitioning / in two 
intervals, one for each box. One gets 

m = I ^ t7-.ov < (3) 

There is a gap centered at to in the profile A(t) and 
A (to) = 2 Aconf /^ (see Fig. 12). In the case of multiple 
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jumps in /, the corresponding gaps merge. To resolve 
single jumps at best, we choose 6t < Tresl'i- In addition, 
H must be as large as possible to get higher profiles A(t) 
where confinement is best detected. Therefore we set 

In order to detect these gaps, one must also make sure 
that their minima are higher than h-Browm the "back- 
ground noise" of pure Brownian trajectories. Indeed, if 
not, the depth of the gap will be reduced, and this will 
corrupt detection. This condition reads 2Acon//5 > 4, or 
Aeon/ > 10. Together with Eq. Q), it can be written in 
terms of the measurable quantities D and Dm'- 

D> 20 Dm- (5) 

In other words, long-term and short-term time scales 
must be well separated. This is observed in a large ma- 
jority of the experimental trajectories below. 

Now we evaluate the capabilities of our algorithm 
on numerical trajectories. They simulate Brownian 
molecules evolving in a square grid of rigid barriers of 
periodicity L. When a step crosses a barrier, it is al- 
lowed with a probability p suitably defined so that the 
long-term diffusion constant equals Dm- Numerical pa- 
rameters match those of experimental trajectories (see 
below). One image is sampled every 40 ms. In addition, 
D is allowed to vary slowly with time in a given trajec- 
tory, to mimic possible composition or physical changes 
of the underlying membrane. More precisely, every sec- 
ond, D is multiplied by a factor randomly chosen in the 
interval [0.9,1.1]. To calculate A, we measure the diffu- 
sion constant Dm by calculating the MSD on intervals 
of 5 s and fitting the slope at the origin (first 3 points). 
We check that Dm — D. Because of statistical fluctua- 
tions on finite samples, all jumps cannot be detected and 
there are false detections. We denote by a the fraction of 
jumps successfully detected by the algorithm among real 
jumps and by ct the ratio of false detections to real jumps. 
The higher a and the lower a, the best the algorithm. To 
localize jumps, we need to estimate the value Aconf since 
gaps are intervals where A(i) is significantly smaller than 
the confinement value Aconf for a sufficient duration. In 
practice, we proceed as follows. We compute Dm as de- 
tailed above. We measure Ar^ and we calculate A(t) 
(see Fig. 121. To avoid biases due to the slow variations of 
D, we calculate the average A{t) of A(t) over successive 
segments of 30 s; A{t) is our estimation of the confined 
profile: A(t) ~ Aconf, because we anticipate that jumps 
are rare. Next we detect intervals where the signal is well 
below A. More precisely, we require that A(t) < aA{t) 
for a duration larger than tc, where the parameters tc 
and a e [0,1] must be optimized. We have scanned large 
ranges of values of both a and tc, and calculated a and 
CT in each case (10'^ trajectories of T =120 s), with the 



typical parameters of the experimental trajectories be- 
low: D ~ 0.1 nm^s~^, L ~ 0.3 /im, Trcs = 10-20 s. 
We observed that tc = St/2 and a — 0.7 gives the best 
compromise with a > 63% and W < 0.7%. The algorithm 
detects two thirds of the jumps and makes very few false 
detections. The value of a comes from the fact that close 
jumps cannot be resolved and are counted only once. 

Another serious complication can arise in experimental 
trajectories. The confinement domains are not necessar- 
ily squares. If they are elongated, like rectangles or more 
complex shapes, Eq. Q is no longer valid. Consider for 
instance rectangles L x pL. There are two types of jumps, 
over short or long edges. It can easily be quantified how 
this affects the relative depths of the gaps associated with 
each type of jump. If the rectangle is extremely elon- 
gated, then only jumps over short edges can be reason- 
ably detected. This complication can be overcome by 
measuring Ax^ and A?/^ and multiplying the latter by 
a counterweight: Acu,r^ = Ax'^ + Ay^/p^. Then both 
kinds of gaps again have the same depth and the pre- 
vious analysis becomes valid. If the main axis of the 
box are not parallel to the axes Ox and Oy, then before 
applying counterweights, one must recover the average 
directions of these main axes by diagonalizing the corre- 
lation matrix C = (r(t)r^(t)) — (r(t))(r'^(t)) averaged 
over sufficiently long time intervals (typically 30 s). We 
checked that this procedure is operational, even though it 
increases significantly the number of false detections be- 
cause of additional numerical operations. However, this 
question goes beyond the scope of this Letter because in 
the experimental trajectories below, by diagonalizing C, 
we find p ~ 1.4 on average (while p ~ 2.1 for a pure un- 
confined random walk [ij]), in which case both types of 
gaps have typically the same depth and a and a are not 
significantly affected if one uses the original profile |^: 
cr > 55% and a < 2.5%. 

Now we apply our algorithm to the 102 experimen- 
tal trajectories from Ref. 0. These are trajectories of 
/i-opioid receptors at the surface of normal rat kidney fi- 
broblast cells, tracked by SPT, after being labeled by 
40 nm gold colloids, at 40 ms time resolution during 
T = 120 s. The parameters D, Dm and L are mea- 
sured II by fitting the MSD by MSD{t) = (LV3)(1 - 
exp(-12Dt/L2)) + ^Dm t. Typically, D ~ 0.1 ixm^s''^, 
L ~ 0.3 pm, Tres ~ 10 s. From Dm and L we deduce Tres- 
is it was noticed that about 15% of the trajectories 
show significant slow variations of D and in a same 
trajectory, up to one order of magnitude for D. This 
could be a serious concern, because variations of these 
parameters cause variations of the reference value Aconf 
that could be misinterpreted as jumps. Fortunately, it 
was also noticed that on individual trajectories, D (x 
when D and L vary along a trajectory. Consequently 
Aconf varies only moderately (see Fig. [3 bottom). This 
is well mimicked by the slow variations of D at fixed 
L that we have imposed in numerical trajectories. If 
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FIG. 2: (Color online) Top: Profile A{t) for a numerical tra- 
jectory with Trea — 20 s. We have also plotted the threshold 
aA{t) with a = 0.7. Real jumps are represented by dia- 
monds and detected ones by circles. One can see 5 detected 
jumps among which one double jump (detected only once). 
Note that the first and last 5t/2 segments are not taken into 
account because A{t) cannot be calculated there. Bottom: 
Upper plot: Profile for an experimental trajectory, with 3 de- 
tected jumps. Here Trea = 9.6 s and the fence model would 
predict 12 jumps. Some intervals where A{t) < aA{t) are 
not considered as jumps because they are not long enough. 
Lower plot: Profile \{t) calculated with Simson et al.'s algo- 
rithm and the threshold Ac = 3.16 (see |l3| for details). 
The receptor is apparently hardly ever confined. 



St — Tres/i exceeds 15s, we set 5t = 15 s, not to loose 
too many points at the beginning and the end of the 
trajectory (see Fig. |2J). We eliminate the 18 trajectories 
that do not satisfy condition 0, as well as those that 
were qualified of "slow or directed diffusion" in Ref. 0, 
because their MSD were more correctly fitted by the cor- 
responding theoretical equations. We are left with 67 
trajectories. 

First of all, we check that A 4 on all profiles. This 
confirms that all trajectories are confined. From the 
value of Tresi we estimate the expected number of jumps, 
if the trajectories were correctly described by a fence or 
corral model, namely (T — 5t) / Tres ■ Then we count the 
detected jumps. An example is provided in Fig. [3 (bot- 
tom). We find that the average ratio of detected jumps 
to the ones expected with fence or corral models is only 
fexp — 16.4%, where we expected more than 55%. The 
histogram in Fig. Ogives greater details. Therefore, as al- 
ready concluded using independent arguments in Ref. , 
a fence or corral model is not able to account alone for 
experimental observations. 

We clearly see that two populations of receptors 
emerge in the histogram. The first population (empty 



bars) contains 19 trajectories, of average detection ratio 
^'exp = 58%, the detected jumps of which can perfectly 
be accounted for by a fence or picket model. The sec- 
ond population (left shaded bar) concerns 48 trajectories 
where we do not detect any jump. Their long-term diffu- 
sion cannot be explained by any fence or picket model. 
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FIG. 3: (Color online) Histogram: ratios of detected jumps 
to the ones expected in a fence model, on experimental 67 
trajectories of /x-opioid receptors An event at a = 400% 
is not shown. 

We are led to the following original conclusion: hop 
diffusion probably exists in trajectories of /i-opioid recep- 
tors, and it can satisfyingly explain the long-term diffu- 
sive behavior of nearly 30 % of the analyzed trajectories. 
But another mechanism must be invoked to explain long- 
term diffusion in a majority of cases. This reinforces the 
need for an alternative model to account for long-term 
diffusion, as proposed in Ref. 0. We might even ask if 
two distinct mechanisms are not independently at work 
in cells to achieve confinement, and if there would not 
exist two populations of /i receptors: the first ones con- 
fined by fences or pickets, and the other ones by long- 
range inter-protein interactions. Additional experiments 
will be necessary to test this hypothesis. 
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